Differences between phytophagous and predatory species in Pentatomidae based on the mitochondrial genome

Abstract Pentatomidae includes many species of significant economic value as plant pests and biological control agents. The feeding habits of Pentatomidae are closely related to their energy metabolism and ecological adaptations. In this study, we sequenced the mitochondrial genomes of 12 Asopinae species using the next‐generation sequencing to explore the effect of dietary changes on mitochondrial genome evolution. Notably, all sequences were double‐stranded circular DNA molecules containing 37 genes and one control region. We then compared and analyzed the mitochondrial genome characteristics of phytophagous and predatory bugs. Notably, no significant difference was observed in the length of the mitochondrial genomes between the predatory and phytophagous bugs. However, the AT content was higher in the mitochondrial genomes of phytophagous bugs than that of predatory bugs. Moreover, phytophagous bugs prefer codon usage patterns ending in A/T compared with predatory bugs. The evolution rate of predatory bugs was lower than that of phytophagous bugs. The phylogenetic relationships across phytophagous bugs' lineages were largely consistent at depth nodes based on different datasets and tree‐reconstructing methods, and strongly supported the monophyly of predatory bugs. Additionally, the estimated divergence times indicated that Pentatomidae explosively radiated in the Early Cretaceous. Subsequently, the subfamily Asopinae and the genus Menida diverged in the Late Cretaceous. Our research results provide data supporting for the evolutionary patterns and classification of Pentatomidae.

A typical insect mitochondrial genome is a circular doublestranded DNA molecules, consisting of 37 genes (13 protein-coding genes (PCGs), 22 transport RNA genes (tRNAs), and two ribosomal RNA genes (rRNAs)), and one control region (Roger et al., 2017;Wang et al., 2015).Currently, insect mitochondrial genomes are widely used for species identification, population genetics, and phylogenetic analysis (Chen, Zheng, et al., 2020;Vico et al., 2020;Wang et al., 2017).Furthermore, we can test traditional classification systems and systematically understand the evolution of classification by analyzing and studying the mitochondrial genomes of different species.
To date, many studies have investigated the phylogenetic relationships of the bugs.Jiang (2017)  Nevertheless, clear and robust evidence exists of Pentatomidae monophyly, involving most of the currently assigned species in the family.Moreover, cyrtophorides are proposed to belong to an independent lineage and be upgraded to Cyrtophoridae.Genevcius et al. (2021), through their study of the tribe Chlorocorini (Pentatominae) using combined DNA and morphological data, revealed that this tribe is not monophyletic.Although phylogenetic studies, including those on representatives of Pentatomidae (Lian et al., 2022;Roca-Cusachs et al., 2022;Zhao, Zhao, et al., 2019), provide a basic framework, phylogenetic relationships within Pentatomidae remain unclear.
In the present study, we sequenced the complete mitochondrial genomes of 12 Asopinae species.We then compared the mitochondrial genomes of phytophagous and predatory bugs, constructed phylogenetic trees, and evaluated the divergence time of Pentatomidae.Our findings could be beneficial for a better understanding of the evolutionary patterns of Pentatomidae and provide a basic theoretical basis for research on biodiversity and biological control.
Specimens were identified based on their morphological characteristics.All specimens were initially stored in 100% ethanol at −20°C, prior to DNA extraction.Voucher specimens were deposited at the Institute of Entomological, Shanxi Agricultural University, Taigu, Shanxi, China.Total genomic DNA was extracted from the thoracic tissue using the Genomic DNA Extraction Kit (Sangon Biotech, Shanghai, China).

| Sequencing, assembly, annotation, and bioinformatics analysis
A whole genome shotgun strategy was employed to construct libraries, that were paired-end (PE 250) sequenced using the Illumina MiSeq sequencing platform.The Fastp software (Chen et al., 2018) was used to obtain high-quality data.The Geneious v.11.0 software (Kearse et al., 2012) was used for sequence assembly and annotation.The reference sequence Arma custos (Fabricius, 1794;NC_051562;Wu et al., 2020) used for the assembly and annotation of the mitogenome of each species was obtained from the NCBI database.The PCGs were identified by open reading frame (ORF; http:// www.ncbi.nlm.nih.gov/ gorf/ gorf.html) using invertebrate mitochondrial genetic codes.The clover secondary structures of the transfer ribonucleic acids (tRNAs) were predicted using the MITOS web server (http:// mitos.bioinf.unileipz ig.de/ ; Bernt et al., 2013).The boundaries of the two rRNAs genes were determined by comparing them with other published rRNA genes in Pentatomidae.The circular maps of the Asopinae mitogenomes were generated using the CGView Server (https:// proks ee.ca/ proje cts/ new).
The PCGs of 64 Pentatomidae species were extracted using Geneious v.11.0, and the amino acid sequences for protein secondary structures were predicted using the SOPMA online website (https:// npsa.lyon.inserm.fr/ cgi-bin/ npsa_ autom at.pl? page=/ NPSA/ npsa_ sopma.html; Geourjon & Deleage, 1995).The codon usage and nucleotide composition of these PCGs were statistically analyzed using MEGA v.11.0 (Tamura et al., 2021).AT and GC skew were calculated as follows: AT (Perna & Kocher, 1995).The effective number of codons (ENC) values, which are commonly used to measure codon bias, of 13 PCGs were calculated using Codon W1.4.2 (Peden, 2000).In order to study the evolutionary patterns between the mitogenomes of phytophagous and predatory bugs, the non-synonymous substitution rate (Ka) and synonymous substitution rate (Ks) of each PCG were calculated using DnaSP v.6.12.03 (Rozas et al., 2017), and the Ka/Ks values were used to determine whether there were natural selection and mutation pressure acting on the protein coding genes.Datamonkey (http:// www.datam onkey.org/ ) was used to perform selective pressure analysis on the PCGs dataset (Murrell et al., 2012(Murrell et al., , 2013)).The tandem repeat sequence of the control region was obtained using the Tandem Repeats Finder server (http:// tandem.bu.edu/ trf/ trf.html; Benson, 1999).

| Phylogenetic analysis
We analyzed phylogenetic relationships among 64 Pentatomidae species with two Scutelleridae species as outgroups (Table 1).
We performed base substitution saturation analysis and sequence composition heterogeneity analysis on the two datasets to determine the feasibility of the phylogeny before constructing a phylogenetic tree.The base substitution saturation index was calculated using DAMBE v.7.0.35 (Xia & Xie, 2001).Heterogeneity analysis was performed using AliGROOVE 1.0.8(Kück et al., 2014).
The phylogenetic trees were generated using the Bayesian Inference (BI) and Maximum likelihood (ML) methods.BI analyses were performed using MrBayes v.3.2.6 (Ronquist et al., 2012), with the GTR + I + G model.The runs were set for 2 × 10 7 generations, with sampling every 1000 generations.The first 25% of generations were removed as burn-in, when the average standard deviation of split frequencies was below 0.01.The ML trees were reconstructed using IQ-TREE v. 2.2.0 (Minh et al., 2020), and the support values for each node were evaluated using the standard bootstrap (BS) algorithm, which was tested 50,000 times.

| Divergence time estimate
Divergence times of Pentatomidae were estimated using the PCGs dataset with a relaxed clock lognormal model in BEAST 1.8.4 (Drummond & Rambaut, 2007).We adopted the GTR + I + G partitioning model and the calibrated Yule model for the prior tree.

| The structure of Asopinae mitochondrial genome
The mitochondrial genome features were comparatively analyzed using 19 Asopinae species (12 newly and 7 previously reported species).All the mitochondrial genomes were double-strand circular DNA molecules (Figure 1), containing 37 genes (13 PCGs, 22 tRNAs, and two rRNAs) and one control region.The arrangement of 37 genes was consistent with the original gene arrangement of Drosophila yakuba Burla, 1954.The general structural characteristics of the mitochondrial genomes are shown in Table S2.The total length of the mitochondrial genome of Asopinae ranged from 15,479 bp (Z.caerulea) to 19,587 bp (Pic.lewisi).In addition, the nucleotide composition of 19 Asopinae species showed a trend of A > T > C > G with a significant AT bias, with the highest (77.14%) and lowest (71.69%)AT content in Z. caerulea and Pic.griseus, respectively.Moreover, all the mitochondrial genomes exhibited a slightly AT-skew (ranging from 0.08 to 0.12, mean = 0.10) and CG-skew (ranging from 0.12 to 0.19, mean = 0.15; Table S3).
The control region of 19 Asopinae species, located between rrnS and trnI (GAT), was the longest non-coding region (661-4651 bp) in the mitochondrial genome (Table S2).C. horvathi and Pic.viridipunctatus exhibited the highest (81.43%) and lowest (67.35%)AT content in the control region, respectively (Table S3).The statistical analysis of tandem repeats in the control region of Asopinae did not reveal any tandem repeat sequences were found in Pin.sanguinipes and Z. caerulea; however, one to six tandem repeat units were found in the other species (Figure 3).

| Genome sizes
We compared and analyzed the mitochondrial genomes lengths (ranging from 14,000 to 20,000 bp) of 19 and 45 species of predatory and phytophagous bugs, respectively (Figure S3).Among the predatory bugs, Z. caerulea and Pic.lewisi exhibited the shortest (15,479 bp) and longest (19,587 bp) mitochondrial genomes, respectively.Among the phytophagous bugs, Graphosoma rubrolineatum (Westwood, 1837) and Nezara viridula (Linnaeus, 1758) exhibited the shortest (15,633 bp) and longest (16,889 bp) mitochondrial genomes, respectively.Notably, differences in genome length may be attributed to non-coding regions among species.In most species of Pentatomidae, the length of the mitochondrial genome ranged from 15,000 to 17,000 bp.However, no significant difference was observed in the mitochondrial genome length between the predatory and phytophagous bugs.
In addition, we compared the lengths of the PCGs between phytophagous and predatory bugs, and found that the length of nad2 in phytophagous bugs was longer than that in predatory bugs (982 ± 13.47 > 959 ± 8.47), indicating significant differences (Figure 4).We then predicted the secondary structures of the PCGs in Pentatomidae.Our findings revealed alpha helix, extended strand, beta turn, and random coil structures (Figures S1).We also compared the mean percentages of these four structures (Figure 5).The percentage of alpha helices in the proteins encoded by atp8, cox3, and nad2 genes in phytophagous bugs was higher than that in predatory bugs.The percentage of extended strands in the proteins encoded by atp8, cox2, and nad2 genes in predatory bugs was higher than that in phytophagous bugs.The percentage of beta turns in the proteins encoded by cytb and nad2 genes in predatory bugs was higher than that in phytophagous bugs.The percentage of random coils in the proteins encoded by nad1 gene in predatory bugs was higher than that in phytophagous bugs.No significant differences were observed in the secondary structures of other proteins.

| Nucleotide composition
The mitochondrial genomes of phytophagous and predatory bugs exhibited a high AT content (Figure S17).Among the predatory and Erthesina fullo (Thunberg, 1783) exhibited the highest (78.94%) and lowest (73.36%)AT content, respectively, with an average AT content of 76.29%.The AT content of the phytophagous bugs was slightly higher than that of the predatory bugs.In addition, the nucleotide composition of the mitochondrial genomes of phytophagous and predatory bugs exhibited AT-skew and CG-skew (Figure 6).Notably, the AT-skew F I G U R E 1 Mitochondrial genomes maps of Asopinae species in this study.
of phytophagous bugs was higher than that of predatory bugs, but there is no significant difference was found in the GC-skew.
In addition, the conserved overlapping regions (trnH/nad4) were found in predatory bugs but not in phytophagous bugs.The longer gene spacers were found between trnS2 and nad1 (20-35 bp) in Pentatomidae, whereas predatory bugs exhibited longer gene spacers between trnM and nad2.

| Codon usage bias
We conducted a comparative analysis of the relative synonymous codon usage (RSCU) between phytophagous and predatory bugs.
We studied the relationships between the ENC and the total codon GC content (GC), first codon GC content (GC1), second codon GC content (GC2), and third codon GC content (GC3) to further investigate codon usage in Pentatomidae (Figure 8).The ENC of the PCGs in phytophagous and predatory bugs exhibited a strong positive correlation with GC and GC3 and a weak positive correlation with GC1 and GC2.

| Evolution rate
By comparing and analyzing the evolutionary rates of phytophagous and predatory bugs, we found that the 13 PCGs of both exhibited Ks > Ka, and Ka/Ks < 1, indicating the evolution of both phytophagous and predatory bugs under purified selection (Figure 9).Among the 13 PCGs in Pentatomidae, the Ka/Ks ratio of atp8 and cox1 was the highest and lowest, respectively.In addition, among the 13 PCGs, it was found that the synonymous substitutions in predatory bugs were higher than those in phytophagous bugs, and the nonsynonymous substitutions of the predatory bugs were lower than those in phytophagous bugs.

F I G U R E 8
Evaluation of codon bias in phytophagous and predatory species in Pentatomidae.

| Phylogenetic analysis
The saturation analysis indicated no saturation in sequences of the two datasets (Iss < Iss.c, and p < .05),and heterogeneity analyses revealed low heterogeneity in sequences (Figures S19 and S20).
Therefore, the two datasets were considered suitable for subsequent phylogenetic analysis.
The phylogenetic trees obtained using the two methods (BI and

| Divergence time estimation
The BEAST analysis indicated that the divergence time of Pentatomidae was 122.78 Mya (95% HPD: 99.15-146.23Mya;

| D ISCUSS I ON AND CON CLUS I ON S
In this study, we sequenced the complete mitochondrial genomes of 12 Asopinae species using second-generation sequencing technology.No gene rearrangements occurred in the mitochondrial genomes of Asopinae, and the sequences were consistent with those of other published Pentatomidae species (Wang et al., 2015;Yuan et al., 2015).The mitochondrial genome size of Pentatomidae is 14-20 kb, with the total lengths in most species ranging from 15 to A comparison of the length of the PCGs and secondary structure of proteins revealed significant differences in the nad2 gene between predatory and phytophagous bugs.Furthermore, phytophagous bugs exhibited a slightly higher AT content than predatory bugs.Moreover, phytophagous bugs tended to prefer codon usage patterns ending in A/T to those of predatory bugs.Different amino acids may cause changes in protein function, thereby affecting organisms and their coevolution with the environment (Hernández- et al., 2008;Liu et al., 2010).In addition, we also found significant differences in the use of Leu between the predatory and phytophagous bugs.Notably, these differences may contribute to the changes that occur in the species to adapt to the environment.The main factors affecting codon bias were mutation pressure and natural selection, with natural selection being the main factor.As insects evolve, the role of natural selection also increases (Behura & Severson, 2013;Nyayanit et al., 2020;Sang, 2019;Wang et al., 2018).The evolutionary rate Ka/Ks < 1 and the selective pressure analysis of Pentatomidae indicated that they are under purified selection.The evolution rate of atp8 was the fastest, whereas that of cox1 was the slowest, which is consistent with the results of the previous studies (Chen, 2022;Ding et al., 2023;Lian et al., 2022).

Montes
Although our results showed that synonymous substitutions in predatory bugs were higher than those in phytophagous bugs, and the non-synonymous substitutions in predatory bugs were lower F I G U R E 1 0 Phylogenetic relationships inferred by the Bayesian Inference (BI) and Maximum likelihood (ML) method based on the protein-coding genes (PCGs) and PRT datasets.Numbers on nodes are the posterior probabilities (PP).
than those in phytophagous bugs, considering that phytophagous bugs have earlier divergence than predatory bugs, it is expected that this relatively young lineage will accumulate more synonymous mutations rather than non-synonymous mutations compared to older phytophagous lineages.Therefore, we may need more factors in the future to explain this result.
We obtained highly consistent topologies of the phylogenetic  the genus Plautia of the Antestiini, which is temporarily placed in this tribe.Notably, we could not determine the relationship between Antestiini and Nezarini, and further research is warranted in this area.
We included three genera-Neojurtina, Pentatoma, and Placosternumof Pentatomini, which is the non-monophyletic poorly defined tribe.
Neojurtina was temporarily identified as a member of Pentatomini (Rider et al., 2018), and Pentatoma formed a monophyletic clade with strong support (1/100/1/100) in our analysis.Therefore, further evidence is required to determine the phylogenetic position of each Pentatomini member.
Phyllocephalinae has a complicated taxonomic history, with the single most diagnostic character being a distinctively short rostrum that does not or only barely surpasses the procoxae.Our analyses, including four representative species of Phyllocephalini, strongly supported its monophyly, which has also been confirmed by Roca-

Cusachs et al. (2022).
The taxon Podopinae is defined as a monophyletic group based on specific morphological characteristics (Rider et al., 2018) Hoplistoderini, Menidini and Asopinae are grouped into one clade.In the previous research (Lian et al., 2022;Roca-Cusachs et al., 2022), it supports the sister-group relationship between Menidini and Asopinae.Moreover, issues in the classification of the tribes in Asopinae remain unresolved.Although various suprageneric names have been proposed, they are not included in the current formal classification (Rider et al., 2018).In addition, our results supported a sister-group relationship between the genus Picromerus and the genus Eocanthecona, whereas enlarged protibia is their distinguishing feature in the traditional classification (Zhao, 2013).The phylogenetic results of this study will provide a good reference for further research on the taxonomic status of Pentatomidae.
Evaluation of the divergence time of Pentatomidae is beneficial for studying its evolutionary history.Notably, the Cretaceous period may have been an important period for the evolution of this group owing to the emergence of warmer and wetter climate conditions globally, as well as the increase in diversity and ecological expansion of angiosperms during this period (Berendse & Scheffer, 2009;Chaboureau et al., 2014;Liu et al., 2019;Yao et al., 2012).During In this study, we conducted a comparative analysis of the mitochondrial genomes of phytophagous and predatory species in the Pentatomidae to explore their evolutionary patterns and understand their evolutionary history, providing data to support research on phylogeny, biodiversity, and biological control.However, the evolutionary information of the mitochondrial genome has not been fully explored, including the genetic information contained in tRNA genes, rRNA genes, and the control regions, as well as the functions, has not been fully explored and requires further in-depth research.
Additionally, the fossil information points of Pentatomidae need further supplementation.Moreover, characterizing, the mitochondrial genome sequences of more species and combining morphological characteristics and molecular evidence is imperative to further explore Pentatomidae evolution.
The two rRNA genes (rrnL and rrnS) were encoded on the Nchain in Asopinae.51.07% of the conserved sites in rrnL were located TA B L E 1 List of sequences used to reconstruct the phylogenetic relationships within Pentatomidae.

F
Organization of the control region in the mitochondrial genomes of Asopinae.The tandem repeats are showed by the green circle with repeat length inside.The orange boxes indicate the length of the sequence of the control region.F I G U R E 4 Sizes of the protein coding genes between phytophagous and predatory bugs.*, ** and *** indicate significant difference between phytophagous and predatory bugs at p < .05,p < .01 and .001,respectively (Nonparametric Tests).All values are mean ± SEM unless otherwise designated.F I G U R E 5 The mean percentages of Alpha helix, Extended strand, Beta turn, and Random coil between phytophagous and predatory bugs.*, ** and *** indicate significant difference between phytophagous and predatory bugs at p < .05,p < .01 and .001,respectively (Nonparametric Tests).All values are mean ± SEM unless otherwise designated.

F
I G U R E 6 AT skew and GC skew of the mitochondrial genomes of Pentatomidae.F I G U R E 7 Use of codons of phytophagous and predatory species in Pentatomidae.| 11 of 18 DING et al.We identified codons under positive selection in the PCGs dataset based on FUBAR and MEME to further analyze the role and direction of selection as the driving force for mitochondrial PCGs evolution.We found pervasive positive/diversifying selection at seven sites and pervasive negative/purifying selection at 3349 sites in FUBAR, with a posterior probability of 0.9.Moreover, we found pervasive positive/diversifying selection at 149 sites in the MEME, with a p-value threshold of .1.
ML) based on the two datasets (PCGs and PRT) demonstrated highly consistent topologies, and most branches exhibited high posterior probability and bootstrap values (Figure 10).The results showed that the phylogenetic relationships of tribes within Pentatominae are relatively chaotic.Notably, Neojurtina typica Distant, 1921 was the earliest diverging lineage within Pentatomidae, and representatives of Nezarini and Antestiini formed sister-groups.Moreover, Caystrini and Halyini also exhibited sister-group relationships.As well as representatives of Strachiini and (Sephelini + Pentatoma) formed sister-group relationships.The subfamily Phyllocephalinae clustered as a monophyletic group and also exhibited a sister-group relationship with the genus Placosternum.Furthermore, G. rubrolineatum and Dybowskyia reticulata (Dallas, 1851) also exhibited a sister-group relationship, and Scotinophara lurida (Burmeister, 1834) and Catacanthus incarnatus (Drury, 1773) were clustered together.The subfamily Asopinae clustered as a monophyletic group and exhibited a sister-group relationship with the tribe Menidini.The relationships within Asopinae were as follows: ((Zicrona + (Troilus + A rma)) + ((Dinorhynchus + Cazira) + (Picromerus + Eocanthecona))).

Figure 11 )
Figure 11), occurring in the Barremian Stage of the Early Cretaceous period in the Mesozoic Era.As the earliest species to separate, the divergence time of N. typica was 97.87 Mya (95% HPD: 73.48-123.12Mya), occurring in the Cenomanian period of the Late Cretaceous in the Mesozoic Era.The divergence time of the subfamily Asopinae and the genus Menida was 67.84 Mya (95% HPD: 50.94-86.16Mya), occurring in the Cretaceous Maastricht period of the Mesozoic Era.The subfamily Asopinae started diverging at 56.20 Mya (95% HPD: 41.78-71.30Mya), occurring in the Cenozoic Paleogene Paleocene Tannitian period.The divergence time of the genus Cazira and D. dybowskyi was 51.51 Mya (95% HPD: 41.78-71.30Mya), occurring in the Eocene Epoch of the Paleogene in the Cenozoic Era.The divergence time of the genus Arma and T. luridus was 33.04 Mya (95% HPD: 23.78-44.03Mya), occurring in the Neogene Oligocene Ruperian period.The divergence time of the genus Picromerus and the genus Eocanthecona was 45.24 Mya (95% HPD: 33.35-57.83Mya), occurring in the Lutai period of the Neogene Eocene in the Cenozoic Era.The divergence time between Pic. lewisi and Pic.bidens was 1.08 Mya (95% HPD: 0.66-1.59Mya), occurring in the Calabrian Stage of the Quaternary Pleistocene in the Cenozoic Era.

F
I G U R E 9 Evolution rate of predatory and phytophagous species in Pentatomidae.| 13 of 18 DING et al. 17 kb.No significant difference was observed in the mitochondrial genome size, which was mainly determined based on the number and length of non-coding regions of phytophagous and predatory bugs belonging to the family Pentatomidae.

F
I G U R E 11 Chronogram with estimated divergence time based on fixed rate calibration among Pentatomidae using BEAST 1.8.4.Horizontal bars represent 95% credibility intervals of time estimates.Numbers on the nodes indicate the mean divergence times.
the evolution of Pentatomidae, three subfamilies (Pentatominae, Phyllocephalinae, and Podopinae) retained phytophagy, and Asopinae shifted to zoophagy.The divergence of the subfamily Asopinae occurred during the Late Cretaceous period of the Mesozoic Era, with subsequent diversification of the most speciose clades in the Cenozoic Era.Although many species exhibit some degree of specialization, none of the Asopinae species are truly host-specific(De Clercq, 2002).Moreover, basic driving factors and evolutionary processes in Pentatomidae are not fully understood and require further research.